function [elements,processInfo] = InitializeSolutionStep...
    (elements,model,processInfo)
ne = model.countElements();
for i=1:ne
    e = model.getElement(i-1);
    f = e.getFace();
    ue = f.getUHatLocal();
    xi = e.getGaussPoint(0);
    w = e.getGaussWeight(0);
    jac = f.jacXAt(xi);
    detJac = det(jac);
    B = zeros(3,6);
    B = e.computeB(f,B,xi);
    C = elements(i).constitutiveMatrix;
    %compute sigma
    if elements(i).isCrack==0
        elements(i).stress = C*B*ue;
    else
        epsilon = B*ue + elements(i).Gc*elements(i).Xi;
        elements(i).stress = C*epsilon;
    end
    %compute renh at element's crack
    if elements(i).isCrack==1
        %compute TGamma and update Kxixi
        lGamma = elements(i).lGamma;
        u_mu = elements(i).Xi;
        t = elements(i).ft+elements(i).H*u_mu(1);
        n = elements(i).crack_normal;
        TGamma = max(t,0)*n;
        if t < 0
            elements(i).Kxixi = w*elements(i).Ge'*...
                C*elements(i).Gc*abs(detJac);
        else
            Ctilde = [elements(i).H*n(1) 0;...
                elements(i).H*n(2) 0];
            elements(i).Kxixi = w*elements(i).Ge'*...
                C*elements(i).Gc*abs(detJac)+...
                lGamma*Ctilde;
        end
        
        %compute renh
        renh = -lGamma*TGamma;
        renh = renh - w*elements(i).Ge'*...
            elements(i).stress*abs(detJac);
        elements(i).renh = renh;
%         renh
    end
    %update ue
    elements(i).ue_old = ue;
end
%elements = calculateOnsetCriteria(elements,model);
end
